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Abstract 

Several methods are constructed for large-scale electronic structure calculations. Test calculations are carried 
out with up to 10 7 atoms. As an application, cleavage process of silicon is investigated by molecular dynamics 
simulation with 10-nm-scale systems. As well as the elementary formation process of the (lll)-(2 x 1) surface, 
we obtain nanoscale defects, that is, step formation and bending of cleavage path into favorite (experimentally 
observed) planes. These results are consistent to experiments. Moreover, the simulation result predicts an explicit 
step structure on the cleaved surface, which shows a bias-dependent STM image. 

Key words: order-iV electronic-structure theory, nanoscale defect, surface process of silicon, bias-dependent STM image. 



1. Introduction 

Quantum mechanical (electronic structure) calcula- 
tion with 10-nm-scale systems is of great importance 
in the present semiconductor technology. Particularly, 
dynamical simulations are highly desirable so as to ex- 
plore industrial processes. These simulations, however, 
are quite difficult for the present standard methodol- 
ogy, such as the Car-Parrinello method owing to 
its heavy computational cost. So as to overcome the 
difficulty, methodology for large systems, with thou- 
sands of atoms or more, has been focused from 1990's. 
130) Its mathematical foundation is the calculation of 
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one-body density matrix, instead of one-electron eigen- 
states. In this paper, we will show the investigation of 
nanoscale defects formed in cleavage process of silicon 
(0), which is an example of our newly-developed prac- 
tical methods of large-scale electronic structure calcu- 
lations, la la La la la llOli These methods are general 
and are applicable to not only single-atom defects but 
also nanoscale defects. 



2. Theory 

The foundation of the large-scale electronic struc- 
ture calculation was established by W. Kohn (2). Any 
physical quantity X can be given by the one-body den- 
sity matrix p as 

{X) = Tv[pX\ = J J drdr'p{r,r')X(r',r), (1) 

where the density matrix is defined, from occupied one- 
electron eigenstates <j>k(r), as 
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Fig. 1. The computational time for bulk Si or Cu as a func- 
tion of the number of atoms (N), up to 11,315,021 atoms 
(Refs. and this work). The time was measured for elec- 

tronic structure calculation with a given atomic structure. The 
calculations were carried out by the conventional eigenstate 
calculation (EIG) and by our methods for large system; (i) 
Krylov-subspace method (KR) with subspace-diagonalization 
procedure, (ii) variational Wannicr-statc method (VR-WS) and 
(iii) perturbative Wannier-state method (PT-WS). For '1CPU' 
computations, we used single Pentium 4™ processor in 2 GHz. 
Parallel computations were carried out by SGI Origin 3800™ 
(for PT-WS method), Origin 2800™ (for VR-WS method) 
and Altix 3700™ (for KR method). 



of our methods, unlike that of the conventional eigen- 
state calculation, is 'order-TV', or linearly proportional 
to the system size (N). Here and hereafter, all the cal- 
culations of Si are carried out using a typical transfer- 
able Hamiltonian. lfl3) For the present calculation of 
Cu, we used the Slater-Koster type Hamiltonian con- 
structed from in the LMTO theory Q. 

As well as the calculation methods for density ma- 
trix, we also constructed a method for Green function 
of large systems, which is based on Krylov subspace. 




Fig. 2. Bending of cleavage path into the favorite (experimen- 
tally observed) planes, (111) and (110) planes. ^4|) The time 
interval At between the snapshots (a) and (b) and that be- 
tween (b) and (c) are At — 1 ps and 2 ps, respectively. 



OCC. 

= J>k><0*|. (2) 
fc 

The most important fact is that, in case that X is a 
short-range operator, the physical quantity (X) is con- 
tributed only by the short-range component of the den- 
sity matrix, even though the density matrix p(r, r') is 
of long range. Owing to this fact, W. Kohn proved that 
the density-functional (ab initio) theory can be con- 
structed from the short-range component of the density 
matrix, instead of eigenstates (Kohn-Sham orbitals). 

Another methodological foundation for large-scale 
calculation is the transferable Hamiltonians in the 
Slater-Koster (tight-binding) form. They are of short 
range and applicable to various circumstances, e.g., 
crystals, defects, liquid and surfaces. Their success has 
been known for decades and can be founded by the ab 
initio theory, lll5h 

For practical large-scale calculations, we have devel- 
oped a set of theories and program codes with obtain- 
ing density matrix. They are founded by generalized 
Wannier state 10; EJ H E3) or Krylov subspace (Refs. 
iflU : Iflt HfT) and this work). These methods are calcu- 
lation methods for short-range component of the den- 
sity matrix with a given Hamiltonian. A bench mark 
is shown in Fig. 1, in which the computational time 



3. Application to cleavage process in silicon 

As the first application of our method, we focused 
on the cleavage process of silicon. Molecular dynamics 
simulations were carried out with 10-nm scale samples 
or with upto 10 5 atoms. (00) 

Since cleavage is a nonequilibrium process with a 
typical velocity scale of propagation velocity, its dy- 
namical mechanism is essential, (d E3) Although 
cleavage process was simulated with electronic struc- 
ture thus far, llisl ; llil) its investigation was quite 
limited, owing to its small system size of 10 2 atoms. 
The present simulation method enables investigation 
beyond the above limitation. 

Numerical accuracy of the present simulation was 
confirmed, not only in elastic constants and surface 
energy, but also in the following quantities that can 
be measured by cleavage experiments; (i) The critical 
stress intensity factor for cleavage K c was evaluated to 
be K c = O.TMPa-ym, which agrees well with experi- 
mental values of K c = 0.65 - 1.24MPayni . iTH) (ii) 
The cleavage propagation velocity was evaluated to be 
v « 2nm/ps=2km/s, which agrees well with a recent 
experimental value of v = 2.3 ± 0.3km/s. <f2CT^> 
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Fig. 3. Elementary surface reconstruction process in silicon 
cleavage with (111) — (2x1) surface, j^) (a)-(c) : Snapshots of 
process, (d)-(c) : Schematic figures of (d) the buckled (2x1) 
structure that appears in the snapshot (b) and (e) the 7r-bondcd 
(Pandcy) (2x1) structure that appears in the snapshot (c). The 
oval indicates the presence of a doubly occupied surface state. 
(c)-(f): Quantum mechanical analysis of a process in which 
the bonding wavefunction {4>i) between A and B sites changes 
into that between A and C sites (A — B A — C); Figure 
(f) shows the interatomic distances (Lab and (Lac- Figure (g) 
shows the occupation weight on the A, B and C atom sites 
(nA- ns i nc) f° r the wavefunction <f>i, A typical transitional 
state, with (Lab ~ (Lac or n B ~ appears at the time 

marked by arrows in (f) and (g). 

3.1. Bending of cleavage path into favorite planes 




Fig. 4. Stable (experimentally observed) cleavage mode on 
(111) -(2x1) surface. 0) (a)-(c):S ucccssivc snapshots with a 
time interval of approximately Ai — 4ps. (d) : Close up of the 
snapshot (c), in which a step structure with a six-mcmbcrcd 
ring appears. 

surface that was proposed first by Pandey J2ll> and is 
currently well established in both experiment and the- 
ory. \'2'A\ Since a surface atom on the ideal (111) surface 
has one dangling-bond electron, the pairing mecha- 
nism between the nearest neighbor dangling-bond elec- 
trons is the motive force of the surface reconstruction. 
Figure 3 shows the actual reconstruction process with 
quantum mechanical analysis of the wavefunction. The 
resultant surface contains a pair of five- and seven- 
membered rings as the unit of the 7r-bonded (2 x 1) 
structure. The 7r-bonding appears between the B and 
D sites in Fig. 3(e), as indicated by an oval. 

3.3. Step structure and bias- dependent STM image 



Experimentally, the cleavage plane of silicon is the 
(111) plane with a metastable (2 x 1) reconstruction or 
the (110) plane (less favorable). (I22T) The appearance 
of these cleaved surfaces cannot be explained from the 
traditional approach with surface energy. As a fruit- 
ful result of the present 10-nm-scale simulation, the 
experimental preference of cleaved surfaces is repro- 
duced by bending of cleavage path; As a typical result, 
Fig. 2 shows that, even if the cleavage initiates arti- 
ficially on another plane ((001) plane), the cleavage 
path is bent into the experimentally observed planes 
((111) and (110) planes). Moreover, we obtained a well- 
defined surface reconstruction on the (111) and (110) 
planes after bending (4). 

3.2. Stable cleavage mode on (111) surface with (2x 1) 
reconstruction 

Our simulation results reproduce the stable (experi- 
mental) cleavage mode with the 7r-bonded (lll)-(2xl) 



In our results, step formation were frequently ob- 
served on the (lll)-(2 x 1) surface. The frequently 
observed step appears in Fig. 4(d) and has a six- 
membered ring at the step edge. The formation pro- 
cess was investigated with quantum mechanical (elec- 
tronic) freedom and we found the mechanism that 
explains why this kind of step appears so frequently. 
(0) Among STM experiments on cleaved surface, on 
the other hand, several step structures were observed 
but their explicit atomic structures have not been set- 
tled. \l'A l24t I25I) Since our results were independent 
from the experiments, we propose the present step 
structure as a hopeful candidate. 

So as to compare with STM experiments, the local 
density of states (LDOS) was calculated for the sur- 
face atoms, the A, B, C, D sites in Fig. 4(d), u sing the 
Krylov subspace method for Green function, (lid : 111!) 
In result, the upper (vacuum side) atom (A or C site) 
has an occupied surface state and the lower (bulk side) 
atom (B or D site) has an unoccupied surface state. The 
present LDOS result corresponds to a bias-dependent 
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STM image. Actually, such a bias-dependent STM im- 
age is observed experimentally J26I) on the flat (non- 
stepped) region, the A and B sites in Fig. 4(d). There- 
fore, we concluded that such a bias-dependent STM 
image will be observed also in the stepped region, the 
C and D sites in Fig. 4(d). More quantitative discus- 
sions will appear elsewhere. 



4. Summary 

Practical methods of large-scale electronic structure 
were constructed and their foundation is to calculate 
(one-body) density matrix or Green function, instead 
of eigenstates. In the application of cleavage process 
of silicon, nanoscale defects, step and bending, were 
observed as well as elementary surface reconstruction 
processes. All the results are consistent to experiments 
qualitatively and quantitatively. Moreover, our simula- 
tion results predict a practical step structure with the 
bias dependency of its STM image. 

Since the present method for large systems is a gen- 
eral quantum mechanical calculation, it has wide ap- 
plications, not specific to silicon or cleavage. Processes 
of other nanoscale or 10-nm-scale systems are possible 
targets in future study. 
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